#include <iostream>
#include <ostream>
#include <istream>
#include <sstream>
#include <fstream>
#include <algorithm>
#include <map>
#include <vector>
#include <deque>
#include <boost/format.hpp>
#include "bcpp_tree.h"
#define _MY_STL_MACROS
#include "stlhelper.h"
#define DEBUG_LEVEL 1
#include "debug_func.h"

using std::string;
using std::cerr;
using std::cout;
using std::vector;
using std::deque;
using std::ostream;
using std::ofstream;
using std::fill;
using std::istream;
using std::stringstream;
using boost::format;


namespace bcpp_tree
{
std::string enode_type_names[] = {	"None",						// eNONE
									"Orphan Tax1 Node",         // eTAX1
									"Orphan Tax2 Node",         // eTAX2
									"Orphan Tax1 Gene",         // eTAX1_GENE
									"Orphan Tax2 Gene",         // eTAX2_GENE
									"Ancestral Node",           // eTAX1_TAX2
									"Ortholog Tax1 Gene",       // eTAX1_GENE_ORTH
									"Ortholog Tax2 Gene",       // eTAX2_GENE_ORTH
									"Ortholog Tax1 Node",       // eTAX1_ORTH
									"Ortholog Tax2 Node",       // eTAX2_ORTH
									"Ortholog Node"};           // eTAX1_TAX2_ORTH


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	t_node

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
t_node::t_node(t_tree* tree_, unsigned i, bool is_leaf, const string& name)
		:	tree(tree_), length(0.0), len_to_root(0.0), len_to_leaf(0.0),
			level(0), node_type(eNONE),
			parent(0), leaf(is_leaf), index(i)
{
	set_name(name);
}

void t_node::print(ostream& os) const
{
	os  << "NODE DUMP"
	<< "\n\tID    = "   << get_name()
	<< "\n\ttype  = "   << enode_type_names[get_node_type()]
	<< "\n\tlevel = "   << get_level()
	<< "\n\tlabel = "   << get_label()
	<< "\n\tlength= "   << get_len_to_root()
	<< "\n";
}

void t_node::set_parent(t_node& new_parent)
{
	// make sure parent is different
	if (&new_parent == parent)
		return;
//	cerr << "Reparent ";
//	if (parent)
//		cerr << "from "	<< parent->get_name();
//	cerr << " to "	<< new_parent.get_name() << "[" << unsigned(this) << "]\n";

	if (parent)
	{
		assert(includes(parent->children, this));
		parent->children.erase(find(beg_to_end(parent->children), this));
	}
	new_parent.children.push_back(this);
	parent = &new_parent;
}

/*
vector<unsigned> t_node::vec_u_id;
t_node* t_node::create_root()
{
	t_node* new_root_node = new t_node(0);
	new_root_node->tree_id = vec_u_id.size();
	vec_u_id.push_back(0);
	++vec_u_id[new_root_node->tree_id];
	return new_root_node;
}

t_node* t_node::add_leaf(const string& name)
{
	t_node* new_node = new t_node(this);
	new_node->set_name(name);
	return new_node;
}
t_node* t_node::add_leaf()
{
	t_node* new_node = new t_node(this);
	return new_node;
}

*/


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	t_tree

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	Operations

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



//________________________________________________________________________________________

//	do_set_len_to_root

//________________________________________________________________________________________
void do_set_len_to_root::operator()(t_node* root)
{
	// set lengths from root
	max_len_to_leaf = 0.0;
	recurse_len_to_root(root);

	// set lengths from leaf
	recurse_len_to_leaves(root);
}

void do_set_len_to_root::recurse_len_to_root(t_node* node)
{
	double len_to_root = node->get_length();
	if (node->get_parent())
		len_to_root += node->get_parent()->get_len_to_root();
	node->set_len_to_root(len_to_root);

	// recurse each child
	if (node->is_leaf())
		// save leaf length
		max_len_to_leaf = std::max(max_len_to_leaf, len_to_root);
	else
		for (unsigned i = 0; i < node->get_children_count(); ++i)
			recurse_len_to_root(node->get_child(i));
}

void do_set_len_to_root::recurse_len_to_leaves(t_node* node)
{
	node->set_len_to_leaf(max_len_to_leaf - node->get_len_to_root());

	// recurse each child
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse_len_to_leaves(node->get_child(i));
}

//________________________________________________________________________________________

//	do_set_node_types

//________________________________________________________________________________________
void do_set_node_types::recurse(t_node* node)
{
	// start from bottom
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));

	// leaves
	if (node->is_leaf())
	{
		if (node->get_name().find(tax1_matching_str) != string::npos)
			node->set_node_type(eTAX1_GENE);
		else if (node->get_name().find(tax2_matching_str) != string::npos)
			node->set_node_type(eTAX2_GENE);
		else
		{
			cerr	<<	"Error:\nUnknown gene:[" <<	node->get_name() << "]\n\n";
			cerr	<< node->get_name()
					<< " , gene 1 : " << tax1_matching_str
					<< " , gene 2 : " << tax2_matching_str
					<< "\n";
		}

		return;
	}


	// get list of leaf node types
	unsigned node_types[eTAX1_TAX2_ORTH + 1];
	fill(node_types,  node_types + eTAX1_TAX2_ORTH + 1,  0);
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		node_types[node->get_child(i)->get_node_type()]++;

	// divide into orthologs, tax1 and tax2 on the basis of leaf nodes

	// = Ancestral_Node
	//		if child includes Ancestral_Node or Ortholog_Node
	if (node_types[eTAX1_TAX2] || node_types[eTAX1_TAX2_ORTH])
		node->set_node_type(eTAX1_TAX2);

	// = Orphan Tax1 Node
	//		if child includes only "Orphan Tax1 Gene" or "Orphan Tax1 Node"
	else if ( (node_types[eTAX1] + node_types[eTAX1_GENE] ) &&
			  !(node_types[eTAX2] + node_types[eTAX2_GENE] ))
		node->set_node_type(eTAX1);

	// = Orphan Tax2 Node
	//		if child includes only "Orphan Tax2 Gene" or "Orphan Tax2 Node"
	else if (!(node_types[eTAX1] + node_types[eTAX1_GENE] ) &&
			 (node_types[eTAX2] + node_types[eTAX2_GENE] ))
		node->set_node_type(eTAX2);

	// = Ortholog_Node
	//		if child includes both tax1 and tax2
	else
	{
		node->set_node_type(eTAX1_TAX2_ORTH);
		tag_ortholog_decendents(node);
	}
}
void do_set_node_types::tag_ortholog_decendents(t_node* node)
{
	// tag orthologous duplication nodes as lineage specific duplications
	// rather than loses
	if (node->get_node_type() == eTAX1)
		node->set_node_type(eTAX1_ORTH);
	else if (node->get_node_type() == eTAX2)
		node->set_node_type(eTAX2_ORTH);
	else if (node->get_node_type() == eTAX1_GENE)
		node->set_node_type(eTAX1_GENE_ORTH);
	else if (node->get_node_type() == eTAX2_GENE)
		node->set_node_type(eTAX2_GENE_ORTH);
	assert (node->get_node_type() != eNONE);

	// recurse
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		tag_ortholog_decendents(node->get_child(i));
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	tree definitions

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



//________________________________________________________________________________________
//	c'tor
//________________________________________________________________________________________
t_tree::t_tree() :
					root(0) ,
					current_node(0),
					current_parent(0),
					finished(false),
					max_level(0)
{
	add_node(false, "root");
	root			=
	current_node    =
	current_parent  = &shared_mem.back();
}


//________________________________________________________________________________________
//	join_nodes
//________________________________________________________________________________________
t_node& t_tree::join_nodes(t_node& node1, t_node& node2)
{
	t_node& new_parent = add_internal_node();
	reparent_to_grandparent();
//	cerr <<
//			format("Merge: [%1$s] %|30t|[%2$s] => %|45t|[%3$s]\n")
//			%	node1.get_name()
//			%	node2.get_name()
//			%	new_parent.get_name();
	node1.set_parent(new_parent);
	node2.set_parent(new_parent);
	return new_parent;
}


//________________________________________________________________________________________
//	add_leaf
//________________________________________________________________________________________
t_node& t_tree::add_node(bool is_leaf, const string& str, unsigned leaf_index)
{
	assert (!finished);

	// do this BEFORE creating node as it might access these parameters
	//	in its constructor
	node_labels.push_back(string());
	node_names.push_back(string());

	shared_mem.push_back(t_node(this, shared_mem.size(), is_leaf, str));

	// save index of leaf
	if (is_leaf)
	{
		leaf_indices[shared_mem.size()-1] = leaf_index;
		leaves.push_back(&shared_mem.back());
	}



	current_node    = &shared_mem.back();
	if (current_parent)
		current_node->set_parent(*current_parent);
	return *current_node;
}


t_node& t_tree::add_leaf(const string& str, unsigned leaf_index)
{
	return add_node(true, str, leaf_index);
}

t_node&  t_tree::add_internal_node()
{
	// internal name
	stringstream os;
	os << "internal_node{" << shared_mem.size() << "}";

	add_node(false, os.str());
	current_parent = current_node;
	return *current_node;
}


//________________________________________________________________________________________
//	finish
//________________________________________________________________________________________
bool t_tree::finish()
{
	debug_func_cout(1);
	finished = true;
	if (current_parent != root)
		return false;

	// make sure levels = parent + 1
	do_set_levels levels;
	for_each(levels);
	max_level = levels.max_level;

	do_sort_by_name sort_by_name;
	for_each(sort_by_name);

	// calculate lengths from root and leaves = parent + 1
	do_set_len_to_root set_len_to_root;
	for_each(set_len_to_root);
	return true;
}

//________________________________________________________________________________________
//	label_current_node
//________________________________________________________________________________________
void t_tree::label_current_node(const string& str)
{
	debug_func_cout(1);
	assert (!finished);
	if (current_node)
		current_node->set_label(str);
	else
		current_parent->set_label(str);
}

//________________________________________________________________________________________
//	reparent_to_grandparent()
//________________________________________________________________________________________
void t_tree::reparent_to_grandparent()
{
	debug_func_cout(1);
	assert (!finished);
	assert(current_parent);
	current_node    = current_parent;
	current_parent  = current_parent->parent;
}

//________________________________________________________________________________________
//	set_branch_len
//________________________________________________________________________________________
void t_tree::set_branch_len(const string& str)
{
	assert (!finished);
	double len = atof(str.c_str());
	current_node->set_length(len);
}

//________________________________________________________________________________________
//	set_branch_len
//________________________________________________________________________________________
void t_tree::set_branch_len(double len)
{
	assert (!finished);
	current_node->set_length(len);
}



//________________________________________________________________________________________
//	reset
//________________________________________________________________________________________
void t_tree::reset()
{
	debug_func_cout(1);
	root			=
	current_node    =
	current_parent  = 0;
	node_names.clear();
	node_labels.clear();
	leaf_indices.clear();
	shared_mem.clear();
	leaves.clear();
	finished		= false;
	add_node(false, "root");
	root			=
	current_node    =
	current_parent  = &shared_mem.back();
}


//________________________________________________________________________________________
//	to_matrix
//________________________________________________________________________________________
void t_tree::to_matrix(t_matrix& matrix)
{
	matrix.allocate(get_leaf_count());
	do_get_matrix get_matrix(matrix);
	for_each(get_matrix);
}





//________________________________________________________________________________________

//	do_newhampshire_format_str

//________________________________________________________________________________________
void do_newhampshire_format_str::operator()(t_node* root)
{
	os << '(';
	for (unsigned i = 0; i < root->get_children_count()-1; ++i)
	{
		recurse(root->get_child(i));
		os << ",";
	}
	recurse(root->get_child(root->get_children_count()-1));
	os << ")\n";
}

void do_newhampshire_format_str::recurse(t_node* node)
{
	unsigned cnt_children = node->get_children_count();
	// is leaf
	if (!cnt_children)
	{
		os << format("%1$s:%2$.5f")
			% node->get_name()
			% node->get_length();
		return;
	}
	os << '(';

	// handle all but last child
	for (unsigned i = 0; i < cnt_children - 1; ++i)
	{
		recurse(node->get_child(i));
		os << ',';
	}

	// handle last child
	recurse(node->get_child(cnt_children-1));

	os << format("):%1$.5f")
			% node->get_length();
}




//________________________________________________________________________________________

//	do_print_tree_full

//________________________________________________________________________________________
void do_print_tree_full::operator()(t_node* root)
{
	recurse(root);
}

void do_print_tree_full::recurse(t_node* node)
{
	os << format("%1$s--[%2$.3f / %3$.3f / %4$.3f]-->%5$s  {%6$s}\n")
		% string((node->get_level()) * 5, ' ')
		% node->get_length()
		% node->get_len_to_root()
		% node->get_len_to_leaf()
		% node->get_name()
		% enode_type_names[node->get_node_type()];
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));
}

//________________________________________________________________________________________

//	do_dump_nodes

//________________________________________________________________________________________
void do_dump_nodes::recurse(t_node* node)
{
	os << format("%1$.3f\t%2$.3f\t%3$s\t%4$s\n")
			% node->get_len_to_leaf()
			% node->get_len_to_root()
			% node->get_name()
			% enode_type_names[node->get_node_type()];
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));
}




//________________________________________________________________________________________

//	do_set_levels

//________________________________________________________________________________________
void do_set_levels::recurse(t_node* node)
{
	for (unsigned i = 0; i < node->get_children_count(); ++i)
	{
		t_node* child = node->get_child(i);
		child->set_level();
		recurse(child);
	}
	max_level = std::max(max_level, node->get_level());
}



//________________________________________________________________________________________

//	do_sort_by_name

//________________________________________________________________________________________
struct op_sort_by_str_ptr
{	bool operator()(const string* a, const string* b) const{return (*a < *b);} };


const string* do_sort_by_name::recurse(t_node* node)
{
	if (node->is_leaf())
		return &node->get_name();

	// only leaves don't have children
	assert(node->get_children_count());

	// get ultimate names of children (name = recursive depth first search to first leaf)
	typedef std::multimap<const string*, unsigned, op_sort_by_str_ptr> t_map_pstr_index;
	t_map_pstr_index names;
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		names.insert(std::pair <const string*, unsigned>(recurse(node->get_child(i)), i));


	// sort children by ultimate name
	std::vector<t_node*> sorted_children;
	for (t_map_pstr_index::iterator i = names.begin(); i != names.end(); ++i)
		sorted_children.push_back(node->children[i->second]);

	assert(stlhelper::is_unique(node->children.begin(), node->children.end()));
	assert(stlhelper::is_unique(sorted_children.begin(), sorted_children.end()));
	sorted_children.swap(node->children);

	return names.begin()->first;

}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	do_collect_orthologs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

//	store_orthologs()


//		store orthologs in the various 1 to many sets

//________________________________________________________________________________________
void do_collect_orthologs::store_genes(t_ortholog_data* store)
{
	//cerr << "{";
	// save total count of genes
	store->cnt_tax1 +=tax1_genes.size();
	store->cnt_tax2 +=tax2_genes.size();

	// save ensg
	store->ensg_sets.push_back(vector<string>());
	store->ensg_length_sets.push_back(vector<double>());
	for (unsigned i = 0; i < tax1_genes.size(); ++i)
	{
		store->ensg_sets.back().push_back(tax1_genes[i]->get_name());
		store->ensg_length_sets.back().push_back(tax1_genes[i]->get_length());
	}
	for (unsigned i = 0; i < tax2_genes.size(); ++i)
	{
		store->ensg_sets.back().push_back(tax2_genes[i]->get_name());
		store->ensg_length_sets.back().push_back(tax2_genes[i]->get_length());
	}

	// clear for next
	tax2_genes.clear();
	tax1_genes.clear();

	//cerr << "}";
}
void do_collect_orthologs::store_orphaned_paralogs()
{
	//cerr << "<";
#ifndef NDEBUG
	// make sure no orthologous genes
	for (unsigned i = 0; i < tax1_genes.size(); ++i)
		assert(tax1_genes[i]->get_node_type() == eTAX1_GENE);
	for (unsigned i = 0; i < tax2_genes.size(); ++i)
		assert(tax2_genes[i]->get_node_type() == eTAX2_GENE);
#endif
	if (!tax1_genes.size() && !tax2_genes.size())
	{
		//cerr << ">";
		return;
	}

	t_ortholog_data* store = 0;


	// tax2 only
	if (!tax1_genes.size())
		store = & orphaned_tax2;

	// tax1 only
	else if (!tax2_genes.size())
		store = & orphaned_tax1;

	// dubious sets
	else
	{
		if (tax1_genes.size() > 1)
		{
			if (tax2_genes.size() > 1)
				store = & orphaned_m_to_m;
			else
				store = & orphaned_m_to_1;
		}
		else
		{
			if (tax2_genes.size() > 1)
				store = & orphaned_1_to_m;
			else
				store = & orphaned_1_to_1;
		}
	}
	store_genes(store);
	//cerr << ">";
}


void do_collect_orthologs::store_orthologs(t_node* node)
{
	//cerr << "(";
#ifndef NDEBUG
	// make sure all orthologous genes
	for (unsigned i = 0; i < tax1_genes.size(); ++i)
	{
		if (tax1_genes[i]->get_node_type() != eTAX1_GENE_ORTH)
		{
			node->print(cerr);
			tax1_genes[i]->print(cerr);
		}
		assert(tax1_genes[i]->get_node_type() == eTAX1_GENE_ORTH);
	}
	for (unsigned i = 0; i < tax2_genes.size(); ++i)
	{
		if (tax2_genes[i]->get_node_type() != eTAX2_GENE_ORTH)
		{
			node->print(cerr);
			tax2_genes[i]->print(cerr);
		}
		assert(tax2_genes[i]->get_node_type() == eTAX2_GENE_ORTH);
	}
#endif
	assert(tax1_genes.size() && tax2_genes.size());

	t_ortholog_data* store = 0;

	if (tax1_genes.size() > 1)
	{
		if (tax2_genes.size() > 1)
			store = & ortho_m_to_m;
		else
			store = & ortho_m_to_1;
	}
	else
	{
		if (tax2_genes.size() > 1)
			store = & ortho_1_to_m;
		else
			store = & ortho_1_to_1;
	}

	store_genes(store);
	//cerr << ")";
}

//________________________________________________________________________________________

//	operator()

//________________________________________________________________________________________
void do_collect_orthologs::recurse(t_node* node)
{
	eNODETYPE node_type = node->get_node_type();
	// take care of adjacent leaves first
	if (node_type == eTAX1_TAX2_ORTH || node_type == eTAX1_TAX2)
		store_orphaned_paralogs();

	// start from the bottom up
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));

	// save list of orthologs in different 1:m, 1:1 etc. lists
	if (node_type == eTAX1_TAX2_ORTH)
		store_orthologs(node);
	else if (node_type == eTAX1_TAX2)
		store_orphaned_paralogs();

	// save actual genes
	else if (node_type == eTAX1_GENE_ORTH || node_type == eTAX1_GENE)
		tax1_genes.push_back(node);
	else if (node_type == eTAX2_GENE_ORTH || node_type == eTAX2_GENE)
		tax2_genes.push_back(node);
}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


//	load_tree_from_newhampshire_format_str


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
boost::tuple<bool, string>
			load_tree_from_newhampshire_format_str(const string& tree_str, t_tree& tree)
{
	string error_str;
	string::const_iterator phylip_iter = tree_str.begin();
	string::const_iterator phylip_end = tree_str.end();
	string str;
	char lastevent = '\0';
	while (phylip_iter != phylip_end)
	{
		switch (*phylip_iter)
		{
			case '(':
				str="";
				tree.add_internal_node();
				lastevent = *phylip_iter;
				break;
			case ')':
			case ',':
				if ( str.length())
				{
					if ( lastevent == ':' )
						tree.set_branch_len(str);
					else
						tree.add_leaf(str, 0);
				}
				else
					tree.add_leaf("", 0);

				if (*phylip_iter == ')')
					tree.reparent_to_grandparent();
				str = "";
				lastevent = *phylip_iter;
				break;
			case ':':
				if (lastevent == '(' || lastevent == ',' )
					tree.add_leaf(str, 0);
				else if (str.length())
					tree.label_current_node(str);
				str = "";
				lastevent = *phylip_iter;
				break;
			case '\n':
			case '\r':
			case '\t':
			case ' ':
				break;
			case ';':
				if (!tree.finish())
				{
					error_str = "Mismatching parentheses";
					return boost::make_tuple(false, error_str);
				}
				return boost::make_tuple(true, error_str);
			default:
				str += *phylip_iter;
		}
		++phylip_iter;
	}
	error_str = "Tree ended prematurely";
	return boost::make_tuple(false, error_str);
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	do_sum_sqrs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


void do_sum_sqrs::operator()(t_node* root)
{
	debug_func_cout(1);
	std::deque<unsigned> parent_branch;
	recurse(root, parent_branch);
	avg_percent_stddev	= cnt_non_null_data == 1 ? 0.0 :
						100.0 * sqrt(sum_sqrs / (cnt_non_null_data * 2.0 - 2.0));
//	APSD =sqrt ( SSQ / (N-2) )    x 100
}



void do_sum_sqrs::recurse(t_node* node, std::deque<unsigned>& parent_branch)
{
	debug_func_cout(1);
	assert(node);

	// save index in parent's list if leaf
	if (node->is_leaf())
	{
		parent_branch.push_back(node->get_leaf_index());
		return;
	}

	assert( node->get_children_count() == 2);

	// create lists of indices for all descendents of each child
	std::deque<std::deque<unsigned> >	branches(node->get_children_count());
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i), branches[i]);

	// length from node to leaf
	double node_dist = node->get_len_to_leaf() * 2;

	assert(branches.size());
	for (unsigned i = 0; i + 1 < branches.size(); ++i)
	{
		for (unsigned j = i + 1; j < branches.size(); ++j)
		{
			std::deque<unsigned>& indices1 = branches[i];
			std::deque<unsigned>& indices2 = branches[j];
			for (unsigned ii = 0; ii < indices1.size(); ++ii)
			{
				for (unsigned jj = 0; jj < indices2.size(); ++jj)
				{
					assert(indices1[ii] < matrix.size());
					assert(indices2[jj] < matrix.size());
					double D = matrix[indices1[ii]][indices2[jj]];
					// ignore null
					if (D <= 0 || D > 5000.0	)
						continue;
					sum_sqrs += 2 * (D - node_dist) * (D - node_dist) / (D *D);
					++cnt_non_null_data;
				}
			}
		}
	}

	// now move all into the branch of my parent
	for (unsigned i = 0; i < branches.size(); ++i)
		parent_branch.insert(parent_branch.begin(), branches[i].begin(), branches[i].end());
	return;
}
void CHECK_LINKS(t_node& curr_node, std::ostream& os)
{
	static unsigned indent = 0;
	os << string(indent, ' ') << curr_node.get_name() << std::endl;
	indent += 2;
	if (!curr_node.is_leaf())
	{
		assert(curr_node.get_child(0));
		assert(curr_node.get_child(1));
		CHECK_LINKS(*curr_node.get_child(0), os);
		CHECK_LINKS(*curr_node.get_child(1), os);
	}
	indent -= 2;
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	op_print_tree

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
void op_print_tree::print() const
{
	tree.for_each(*this);
}

void op_print_tree::print(const t_node& curr_node) const
{
//	debug_func_cout(1);
	os.precision(precision);
	os << std::fixed;
	recurse(curr_node);

}
void op_print_tree::recurse(const t_node& curr_node) const
{
//	debug_func_cout(1);
	if (curr_node.is_leaf())
	{
		print_horizontal_lines();
		os << curr_node.get_name() << std::endl;
		return;
	}
	else
	{
		print_horizontal_lines();
		levels.push_back(1);
		recurse(*curr_node.get_child(0));
		print_vertical_lines();
		os << curr_node.get_len_to_leaf();
		levels.back() = 0;
		os << std::endl;
		print_vertical_lines();
		recurse(*curr_node.get_child(1));
		levels.pop_back();
	}
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	do_get_matrix

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


void do_get_matrix::operator()(t_node* root)
{
	branch_lengths.assign(root->get_tree()->get_leaf_count(), 0.0);

	std::deque<unsigned> parent_branch;
	recurse(root, parent_branch);
}



void do_get_matrix::recurse(t_node* node, std::deque<unsigned>& parent_branch)
{
	assert(node);

	// save index in parent's list if leaf
	if (node->is_leaf())
	{
		assert(branch_lengths.size() > node->get_leaf_index());
		parent_branch.push_back(node->get_leaf_index());
		branch_lengths[node->get_leaf_index()] += node->get_length();
		return;
	}

	assert( node->get_children_count() == 2);

	// create lists of indices for all descendents of each child
	std::deque<std::deque<unsigned> >	branches(node->get_children_count());
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i), branches[i]);

	assert(branches.size());
	for (unsigned i = 0; i + 1 < branches.size(); ++i)
	{
		for (unsigned j = i + 1; j < branches.size(); ++j)
		{
			std::deque<unsigned>& indices1 = branches[i];
			std::deque<unsigned>& indices2 = branches[j];
			for (unsigned ii = 0; ii < indices1.size(); ++ii)
			{
				for (unsigned jj = 0; jj < indices2.size(); ++jj)
				{
					unsigned gene_1 = indices1[ii];
					unsigned gene_2 = indices2[jj];
					assert(gene_1 < matrix.size());
					assert(gene_2 < matrix.size());
					matrix[gene_1][gene_2] =
					matrix[gene_2][gene_1] =
						branch_lengths[gene_2] +
						branch_lengths[gene_1];
				}
			}
		}
	}

	// now move all into the branch of my parent
	for (unsigned i = 0; i < branches.size(); ++i)
		parent_branch.insert(parent_branch.begin(), branches[i].begin(), branches[i].end());
	for (unsigned i = 0; i < parent_branch.size(); ++i)
		branch_lengths[parent_branch[i]] += node->get_length();
	return;
}
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	do_sum_sqrs2

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


void do_sum_sqrs2::operator()(t_node* root)
{
	debug_func_cout(1);
	branch_lengths.assign(root->get_tree()->get_leaf_count(), 0.0);
	branch_differences.assign(root->get_tree()->get_leaf_count(), 0.0);
	std::deque<unsigned> parent_branch;
	recurse(root, parent_branch);
	avg_percent_stddev	= cnt_non_null_data == 1 ? 0.0 :
						100.0 * sqrt(sum_sqrs / (cnt_non_null_data * 2.0 - 2.0));
//	APSD =sqrt ( SSQ / (N-2) )    x 100
}



void do_sum_sqrs2::recurse(t_node* node, std::deque<unsigned>& parent_branch)
{
//	debug_func_cout(1);
	assert(node);

	// save index in parent's list if leaf
	if (node->is_leaf())
	{
		assert(branch_lengths.size() > node->get_leaf_index());
		parent_branch.push_back(node->get_leaf_index());
		branch_lengths[node->get_leaf_index()] += node->get_length();
		return;
	}

	assert( node->get_children_count() == 2);

	// create lists of indices for all descendents of each child
	std::deque<std::deque<unsigned> >	branches(node->get_children_count());

	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i), branches[i]);
	double node_dist = node->get_len_to_leaf() * 2;

	assert(branches.size());
	for (unsigned i = 0; i + 1 < branches.size(); ++i)
	{
		for (unsigned j = i + 1; j < branches.size(); ++j)
		{
			std::deque<unsigned>& indices1 = branches[i];
			std::deque<unsigned>& indices2 = branches[j];
			for (unsigned ii = 0; ii < indices1.size(); ++ii)
			{
				for (unsigned jj = 0; jj < indices2.size(); ++jj)
				{
					unsigned gene_1 = indices1[ii];
					unsigned gene_2 = indices2[jj];
					assert(fabs(node_dist - (branch_lengths[gene_2] +
											 branch_lengths[gene_1])) < 0.001);;

					assert(gene_1 < matrix.size());
					assert(gene_2 < matrix.size());
					double D = matrix[gene_1][gene_2];
					// ignore null
					if (D <= 0.0 || D > 5000.0	)
						continue;
					double diff = 2 * (D - node_dist) * (D - node_dist) / (D *D);
					sum_sqrs += diff;
					branch_differences[gene_1] += diff;
					branch_differences[gene_2] += diff;
					++cnt_non_null_data;

				}
			}
		}
	}

	// now move all into the branch of my parent
	for (unsigned i = 0; i < branches.size(); ++i)
		parent_branch.insert(parent_branch.begin(), branches[i].begin(), branches[i].end());
	for (unsigned i = 0; i < parent_branch.size(); ++i)
		branch_lengths[parent_branch[i]] += node->get_length();
	return;
}
}// namespace



